 // <h1>Modification of Example 4 of libMesh source</h1>
 //
 // This is the fourth example program.  It builds on
 // the third example program by showing how to formulate
 // the code in a dimension-independent way.  Very minor
 // changes to the example will allow the problem to be
 // solved in one, two or three dimensions.
 //
 // This example will also introduce the PerfLog class
 // as a way to monitor your code's performance.  We will
 // use it to instrument the matrix assembly code and look
 // for bottlenecks where we should focus optimization efforts.
 //
 // This example also shows how to extend example 3 to run in
 // parallel.  Notice how litte has changed!  The significant
 // differences are marked with "PARALLEL CHANGE".


// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <math.h>

// Basic include file needed for the mesh functionality.
#include "mesh_generation.h"
#include "exodusII_io.h"
#include "gnuplot_io.h"

#include "linear_implicit_system.h"
#include "projection_errors.h"

// Define the Finite Element object.
#include "fe.h"

// Define the DofMap, which handles degree of freedom
// indexing.
#include "dof_map.h"
#include "quadrature_gauss.h"
#include "parallel.h"

#include "libmesh_project_soln.h"
#include "moab_project_soln.h"

// partitioners
#include "linear_partitioner.h"
#include "metis_partitioner.h"
#include "parmetis_partitioner.h"
#include "sfc_partitioner.h"

// Define the PerfLog, a performance logging utility.
// It is useful for timing events in a code and giving
// you an idea where bottlenecks lie.
#include "perf_log.h"

// The definition of a geometric element
#include "elem.h"

#include "string_to_enum.h"
#include "getpot.h"

// Bring in everything from the libMesh namespace
using namespace libMesh;

// Begin the main program.
int main (int argc, char** argv)
{
	// Initialize libMesh and any dependent libaries, like in example 2.
	LibMeshInit init (argc, argv);

	// Declare a performance log for the main program
	// PerfLog perf_main("Main Program");

	// Create a GetPot object to parse the command line
	GetPot command_line (argc, argv);

	// Check for proper calling arguments.
	if (argc < 2)
	{
		if (libMesh::processor_id() == 0)
			std::cerr << "Usage:\n"
				  <<"\t " << argv[0] << " -d 2(3) -o FIRST(SECOND)" << " -n 15"
				  << std::endl;

		// This handy function will print the file name, line number,
		// and then abort.  Currrently the library does not use C++
		// exception handling.
		libmesh_error();
	}
	else
	{
		// Brief message to the user regarding the program name
		// and command line arguments.
		std::cout << "Running " << argv[0];

		for (int i=1; i<argc; i++)
		std::cout << " " << argv[i];

		std::cout << std::endl << std::endl;
	}


	// Read problem dimension from command line.  Use int
	// instead of unsigned since the GetPot overload is ambiguous
	// otherwise.
	int dim = 3;
	if ( command_line.search(1, "-d") )
		dim = command_line.next(dim);

	// Skip higher-dimensional examples on a lower-dimensional libMesh build
	libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");

	// Create a mesh with user-defined dimension.
	// Read number of elements from command line
	int ps = 2;
	if ( command_line.search(1, "-n") )
		ps = command_line.next(ps);

	bool compute_error = false;
	if ( command_line.search(1, "-error") )
		compute_error = true;

	bool io = false;
	if ( command_line.search(1, "-io") )
		io = true;

	bool disable_moab_projection = false;
	if ( command_line.search(1, "-no_moab") )
		disable_moab_projection = true;

	bool disable_libmesh_projection = false;
	if ( command_line.search(1, "-no_libmesh") )
		disable_libmesh_projection = true;

	bool new_moab_method = false;
	if ( command_line.search(1, "-new_moab_proj") )
		new_moab_method = true;

	// Read FE order from command line
	std::string order = "FIRST";
	if ( command_line.search(1, "-o") )
		order = command_line.next(order);

	bool verbose = false;
	if ( command_line.search(1, "-g") )
		verbose = true;

	// Set FE Family
	std::string family = "LAGRANGE";

	Mesh src_mesh, target_lmesh, target_moab;

	// generate the mesh as requested
	generate_mesh(dim, ps, order, src_mesh);

	// Create an equation systems object.
	EquationSystems src_equation_systems (src_mesh);
	EquationSystems *target_equation_systems_lmesh = NULL;
	EquationSystems *target_equation_systems_moab = NULL;

	// Declare the system and its variables.
	// Create a system named PRIMARY_SYSTEM_NAME
	LinearImplicitSystem& src_system = src_equation_systems.add_system<LinearImplicitSystem> (PRIMARY_SYSTEM_NAME);
	ExplicitSystem& src_system_sub = src_equation_systems.add_system<ExplicitSystem> (AUXILLARY_SYSTEM_NAME);

	// Add the variable "SOLUTION" and "PARTITION" to the system.
	src_system.add_variable("SOLUTION",
					  Utility::string_to_enum<Order>   (order),
					  Utility::string_to_enum<FEFamily>(family));
	src_system_sub.add_variable("PARTITION", CONSTANT, MONOMIAL);

	// Initialize the data structures for the equation system.
	src_equation_systems.init();
	fill_auxillary_systems(src_equation_systems, AUXILLARY_SYSTEM_NAME) ;
	// Print information about the system to the screen.
	src_equation_systems.print_info();
	src_mesh.print_info();

	if (!disable_libmesh_projection)
	{
		// generate the new target mesh
		generate_mesh(dim, ps*2, order, target_lmesh);

		// create the libMesh equation system
		target_equation_systems_lmesh = new EquationSystems(target_lmesh);

		LinearImplicitSystem& target_system_lmesh = target_equation_systems_lmesh->add_system<LinearImplicitSystem> (PRIMARY_SYSTEM_NAME);
		ExplicitSystem& target_system_lmesh_sub = target_equation_systems_lmesh->add_system<ExplicitSystem> (AUXILLARY_SYSTEM_NAME);

		target_system_lmesh.add_variable("SOLUTION",
						  Utility::string_to_enum<Order>   (order),
						  Utility::string_to_enum<FEFamily>(family));
		target_system_lmesh_sub.add_variable("PARTITION", CONSTANT, MONOMIAL);

		target_equation_systems_lmesh->init();

		fill_auxillary_systems(*target_equation_systems_lmesh, AUXILLARY_SYSTEM_NAME) ;

		target_equation_systems_lmesh->print_info();
		target_lmesh.print_info();
	}


	if (!disable_moab_projection)
	{
		// generate the new target mesh
		generate_mesh(dim, ps*2, order, target_moab);

		// create the libMesh-MOAB equation system
		target_equation_systems_moab = new EquationSystems(target_moab);

		LinearImplicitSystem& target_system_moab = target_equation_systems_moab->add_system<LinearImplicitSystem> (PRIMARY_SYSTEM_NAME);
		ExplicitSystem& target_system_moab_sub = target_equation_systems_moab->add_system<ExplicitSystem> (AUXILLARY_SYSTEM_NAME);

		target_system_moab.add_variable("SOLUTION",
						  Utility::string_to_enum<Order>   (order),
						  Utility::string_to_enum<FEFamily>(family));
		target_system_moab_sub.add_variable("PARTITION", CONSTANT, MONOMIAL);

		target_equation_systems_moab->init();

		fill_auxillary_systems(*target_equation_systems_moab, AUXILLARY_SYSTEM_NAME) ;

		target_equation_systems_moab->print_info();
		target_moab.print_info();
	}

	// Perform a standard projection based on analytical function
	// Optionally, we can solve a real system and use that too :)
	//src_equation_systems.get_system(PRIMARY_SYSTEM_NAME).solve();
	ExactPoisson exact_fn;
	src_equation_systems.get_system(PRIMARY_SYSTEM_NAME).project_solution(&exact_fn) ;

	// project the source solution on to the libMesh refined mesh
	if (!disable_libmesh_projection)
	{
    LibMeshProjection libmesh_project(verbose);
    libmesh_project.set_source_equation_system(src_equation_systems);
    libmesh_project.set_target_equation_system(*target_equation_systems_lmesh);

    libmesh_project.init() ;

		libmesh_project.apply () ;
	}

	if (!disable_moab_projection)
	{
    MOABProjection moab_project(verbose);
		moab_project.use_new_moab_method(new_moab_method) ;
    moab_project.set_source_equation_system(src_equation_systems);
    moab_project.set_target_equation_system(*target_equation_systems_moab);

    moab_project.init() ;

		// project the source solution on to the MOAB refined mesh
		moab_project.apply ((io ? "test.h5m" : "")) ;
	}

	if (compute_error)
	{
		// now compute the errors in the solution via our custom functions
		std::vector<Real> src_errors, libMesh_errors, MOAB_errors;

		compute_errors(src_equation_systems, src_errors, PRIMARY_SYSTEM_NAME);

		if (!disable_libmesh_projection)
			compute_errors(*target_equation_systems_lmesh, libMesh_errors, PRIMARY_SYSTEM_NAME);

		if (!disable_moab_projection)
			compute_errors(*target_equation_systems_moab, MOAB_errors, PRIMARY_SYSTEM_NAME);

		std::cout << "\n\n\t\t========================================================================\n" ;
		std::cout << "\t\t|                        Errors in solutions                            |\n" ;
		std::cout << "\t\t========================================================================\n\n" ;
		std::cout << "\t\t|\t\t[source solution]        =   " << src_errors[0] << "\t\t\t|" << std::endl ;
		if (!disable_libmesh_projection)
			std::cout << "\t\t|\t\t[projected via libMesh]  =   " << libMesh_errors[0] << "\t\t\t|" << std::endl ;
		if (!disable_moab_projection)
			std::cout << "\t\t|\t\t[projected via MOAB]     =   " << MOAB_errors[0] << "\t\t\t|" << std::endl ;
		std::cout << "\n\t\t========================================================================" << std::endl << std::endl ;
	}

	// After solving the system write the solution
	if (io)
	{
		if(dim == 1)
		{
			GnuPlotIO srcplot(src_mesh,"Example 4, 1D",GnuPlotIO::GRID_ON);
			srcplot.write_equation_systems("src_mesh_1",src_equation_systems);
			if (!disable_libmesh_projection)
			{
				GnuPlotIO tarplot(target_lmesh,"Example 4, 1D",GnuPlotIO::GRID_ON);
				tarplot.write_equation_systems("target_libmesh_1", *target_equation_systems_lmesh);
			}
			if (!disable_moab_projection)
			{
				GnuPlotIO tarplot_moab(target_moab,"Example 4, 1D",GnuPlotIO::GRID_ON);
				tarplot_moab.write_equation_systems("target_mesh_1", *target_equation_systems_moab);
			}
		}
		else
		{
			ExodusII_IO (src_mesh).write_discontinuous_exodusII ((dim == 3) ?
			"src_mesh_3.e" : "src_mesh_2.e",src_equation_systems);
			if (!disable_libmesh_projection)
				ExodusII_IO (target_lmesh).write_discontinuous_exodusII ((dim == 3) ?
					"target_libmesh_3.e" : "target_libmesh_2.e", *target_equation_systems_lmesh);
			if (!disable_moab_projection)
				ExodusII_IO (target_moab).write_discontinuous_exodusII ((dim == 3) ?
					"target_moab_3.e" : "target_moab_2.e", *target_equation_systems_moab);
		}
	}

	if (!disable_libmesh_projection)
		delete target_equation_systems_lmesh ;

	if (!disable_moab_projection)
		delete target_equation_systems_moab ;

	// All done.
	return 0;
}

void compute_errors(EquationSystems& es, std::vector<Real>& errors, std::string primary_sys_name)
{
  const bool compute_L2_error = true;
  const bool compute_H1_error = false;

  errors.resize(2);

	// Get a constant reference to the mesh object.
	const MeshBase& mesh = es.get_mesh();

	// The dimension that we are running
	const unsigned int dim = mesh.mesh_dimension();

	// Get a reference to the LinearImplicitSystem we are solving
	LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>(primary_sys_name);

	// A reference to the \p DofMap object for this system.  The \p DofMap
	// object handles the index translation from node and element numbers
	// to degree of freedom numbers.  We will talk more about the \p DofMap
	// in future examples.
	const DofMap& dof_map = system.get_dof_map();

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = dof_map.variable_type(0);

  // Build a Finite Element object of the specified type.  Since the
  // \p FEBase::build() member dynamically creates memory we will
  // store the object as an \p AutoPtr<FEBase>.  This can be thought
  // of as a pointer that will clean up after itself.
  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));

  // A 5th order Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, SIXTH);

  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule (&qrule);

  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  // We begin with the element Jacobian * quadrature weight at each
  // integration point.
  const std::vector<Real>& JxW = fe->get_JxW();

  // The physical XY locations of the quadrature points on the element.
  // These might be useful for evaluating spatially varying material
  // properties at the quadrature points.
  const std::vector<Point>& q_point = fe->get_xyz();

  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real> >& phi = fe->get_phi();

  // The element shape function gradients evaluated at the quadrature
  // points.
  const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();

	// This vector will hold the degree of freedom indices for
	// the element.  These define where in the global system
	// the element degrees of freedom get mapped.
	std::vector<unsigned int> dof_indices;

	system.update() ;

	ExactPoisson exact;

	// Now we will loop over all the elements in the mesh.
	// We will compute the element matrix and right-hand-side
	// contribution.  See example 3 for a discussion of the
	// element iterators.
	MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
	const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

	for ( ; el != end_el; ++el)
	{
	  const Elem* elem = *el;

	  // Get the degree of freedom indices
	  dof_map.dof_indices (elem, dof_indices, 0);

    // Compute the element-specific data for the current
    // element.  This involves computing the location of the
    // quadrature points (q_point) and the shape functions
    // (phi, dphi) for the current element.
    fe->reinit (elem);

    for (unsigned int qp=0; qp<qrule.n_points(); qp++)
    {
      Real solution = 0.0, deriv = 0.0;
      for (unsigned int i=0; i<phi.size(); i++)
      {
          solution += phi[i][qp]*(*system.current_local_solution)(dof_indices[i]) ;
          // deriv += dphi[j][qp]);
      }

      Real exact_sol = exact(q_point[qp]) ;

      //std::cout << "\nQP="<< qp << "\tExact="<< exact_sol << "\tComputed="<< solution;

      errors[0] += JxW[qp]*(solution-exact_sol)*(solution-exact_sol);
    }
	}

  // Do an all-reduce to aggregate the error computed
	Parallel::sum(errors) ;

	errors[0] = std::sqrt(errors[0]);
}

void fill_auxillary_systems(EquationSystems& es, std::string aux_sys_name)
{
	// Get a reference to the ExplicitSystem we are updating
	ExplicitSystem& aux_system = es.get_system<ExplicitSystem>(aux_sys_name);

	// Get a constant reference to the mesh object.
	const MeshBase& mesh = es.get_mesh();

	// This vector will hold the degree of freedom indices for
	// the element.  These define where in the global system
	// the element degrees of freedom get mapped.
	std::vector<unsigned int> dof_indices;
	DenseVector<Number> Fe;

  PetscInt rank;
	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

/*
	// Now we will loop over all the elements in the mesh.
	// We will compute the element matrix and right-hand-side
	// contribution.  See example 3 for a discussion of the
	// element iterators.
	MeshBase::const_node_iterator       nd     = mesh.local_nodes_begin();
	const MeshBase::const_node_iterator end_nd = mesh.local_nodes_end();

	std::vector<Real> dofs(mesh.n_local_nodes());
	std::vector<unsigned int> indices(mesh.n_local_nodes());

  unsigned int index = 0;
	for ( ; nd != end_nd; ++nd, ++index)
	{
	  const Node* node = *nd;

	  indices[index] = node->dof_number(0, 0, 0) ;
	  dofs[index] = static_cast<double>(index);

//	  indices[index] = node->dof_number(0, 0, 0) ;
//	  dofs[index] = static_cast<double>(node->id());
	}
*/

  MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();

	std::vector<Real> subdomains(mesh.n_local_nodes());
	std::vector<unsigned int> indices(mesh.n_local_nodes());

  unsigned int index = 0;
	for ( ; el != end_el; ++el, ++index)
	{
	  const Elem* elem = *el;

	  indices[index] = elem->id();
	  subdomains[index] = rank;
//	  PetscPrintf("\n[%D] Element %d", rank, elem->id());
	}

  // will overwrite the shared dofs again but we dont care.
  aux_system.solution->insert (subdomains, indices);

	aux_system.solution->close() ;
	aux_system.update	()	;
//	aux_system.solution->print() ;
//	aux_system.solution->print_global() ;
//	std::cin.get();
}

void generate_mesh(unsigned int dim, unsigned int ps, const std::string& order, Mesh& mesh)
{
  // Use the MeshTools::Generation mesh generator to create a uniform
  // grid on the square [-1,1]^D.  We instruct the mesh generator
  // to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27
  // elements in 3D.  Building these higher-order elements allows
  // us to use higher-order approximation, as in example 3.

  const Real width = dim > 1 ? 1. : 0.;
  const Real height = dim > 2 ? 1. : 0.;

  if (order == "FIRST")
    {
      // No reason to use high-order geometric elements if we are
      // solving with low-order finite elements.
      MeshTools::Generation::build_cube (mesh,
                                         ps,
										 (dim>1) ? ps : 0,
										 (dim>2) ? ps : 0,
                                         0, 1.,
                                         0, width,
                                         0, height,
                                         (dim==1)    ? EDGE2 :
                                         ((dim == 2) ? QUAD4 : HEX8));
    }

  else
    {
      MeshTools::Generation::build_cube (mesh,
										 ps,
										 (dim>1) ? ps : 0,
										 (dim>2) ? ps : 0,
                                         0, 1.,
                                         0, width,
                                         0, height,
                                         (dim==1)    ? EDGE3 :
                                         ((dim == 2) ? QUAD9 : HEX27));
    }

    // repartition the mesh
    SFCPartitioner partn ;
    partn.repartition(mesh);
}


